Doruk Efe Gökmen -- 01/08/2018 -- Ankara
Consider 1D gaussian integral $I$. We calculate it as follows: $I^2=\int^\infty_{-\infty}\frac{\text{d}y}{\sqrt{2\pi}} \int^\infty_{-\infty}\frac{\text{d}x}{\sqrt{2\pi}}e^{-(x^2+y^2)/2}=\int^{2\pi}_0 \frac {\text{d}\phi}{2\pi}\int^\infty_0 r e^{-r^2/2} \text{d}r = \int^{2\pi}_0 \frac {\text{d}\phi}{2\pi}\int^1_0 \text{d}\Upsilon=1$.
The variable transformations indicate a transformation between uniformly distributed random variables $\Upsilon, \phi$ and gaussian distributed random variables $x, y$. Changes of variables in integrals also apply to samples. This is a relation between the integration variables on the exponent of the gaussian and the gaussian distributed random variables.
%pylab inline
import random, math, pylab
def gauss_test(sigma):
#the 2nd transformed random variables are uniformly distributed
phi = random.uniform(0.0, 2.0 * math.pi)
Upsilon = random.uniform(0.0, 1.0)
#the first transformed variables are
Psi = - math.log(Upsilon)
r = sigma * math.sqrt(2.0 * Psi)
#the original variables
x = r * math.cos(phi)
y = r * math.sin(phi)
return [x, y]
#exact Gaussian distrubution:
list_x = [i * 0.1 for i in xrange(-40, 40)]
list_y = [math.exp(- x ** 2 / 2.0) / (math.sqrt(2.0 * math.pi)) for x in list_x]
#sampled distribution:
n_sampled_pairs = 50000
data = []
for sample in xrange(n_sampled_pairs):
data += gauss_test(1.0)
#graphics output
pylab.plot(list_x, list_y, color='b', label='exact')
pylab.hist(data, bins=150, normed=True, color='r', histtype='step', label='sampled')
pylab.legend()
pylab.title('Sampling of the gaussian distribution')
pylab.xlabel('$x$', fontsize=14)
pylab.ylabel('$\pi(x)$', fontsize=14)
pylab.savefig('plot-gauss_test.png')
Populating the interactive namespace from numpy and matplotlib
The plots of independently distributed Gaussian random variables is given below. In the figure, each point is an independent random variable in $x$ and $y$. The resulting 2D distribution is isotropic. The reason for this can be seen in the previous section where we have used the variable transformation to get an integral over angle $\phi$ on the $xy$ plane which can be seen as a randomly distributed number between $0$ and $2\pi$ in the evaluation of the gaussian integral $I$. Hence, this property is unique of gaussian distributions, and is generally valid for gaussians in any dimensions, as the 3 variable case (plotted in the second figure) also suggests so.
%pylab inline
import random, math, pylab, mpl_toolkits.mplot3d
nsamples = 10000
x_list, y_list, z_list, x_list_n2, y_list_n2, x_list_n3, y_list_n3, z_list_n3 = [], [], [], [], [], [], [], []
for sample in range(nsamples):
#two gaussian distributed random variables:
x_list.append(random.gauss(0.0, 1.0))
y_list.append(random.gauss(0.0, 1.0))
z_list.append(random.gauss(0.0, 1.0))
#begin graphics output
#2d distribution
pylab.plot(x_list, y_list, color='black', marker='.', linestyle='')
pylab.title('Samples from the 2D Gaussian distribution')
pylab.xlabel('$x$', fontsize=14)
pylab.ylabel('$y$', fontsize=14)
pylab.xlim(-4.0, 4.0)
pylab.ylim(-4.0, 4.0)
pylab.axes().set_aspect('equal') # set the aspect ratio of the plot
pylab.savefig('plot-gauss_2d.png')
#3d distribution
fig = pylab.figure()
ax = fig.gca(projection='3d')
ax.set_aspect('equal') # set the aspect ratio of the plot
pylab.plot(x_list, y_list, z_list, '.')
pylab.title('Samples from the 3D gaussian distribution')
ax.set_xlabel('$x$', fontsize=14)
ax.set_ylabel('$y$', fontsize=14)
ax.set_zlabel('$z$', fontsize=14)
pylab.show()
Populating the interactive namespace from numpy and matplotlib
%pylab inline
import random, math, pylab, mpl_toolkits.mplot3d
#we can also give the distribution a "haircut"! Let us normalise each point to a radius equal to 1:
nsamples = 5000
x_list_n2, y_list_n2, x_list_n3, y_list_n3, z_list_n3 = [], [], [], [], []
for sample in xrange(nsamples/5):
x, y, z = random.gauss(0.0, 1.0), random.gauss(0.0, 1.0), random.gauss(0.0, 1.0)
radius_2 = math.sqrt(x ** 2 + y ** 2)
x_list_n2.append(x / radius_2)
y_list_n2.append(y / radius_2)
radius_3 = math.sqrt(x ** 2 + y ** 2 + z ** 2)
x_list_n3.append(x / radius_3)
y_list_n3.append(y / radius_3)
z_list_n3.append(z / radius_3)
#figure output
pylab.plot(x_list_n2, y_list_n2, color='black', marker='.', linestyle='')
pylab.title('Uniform sampling of the circle')
pylab.xlabel('$x$', fontsize=14)
pylab.ylabel('$y$', fontsize=14)
pylab.xlim(-2.0, 2.0)
pylab.ylim(-2.0, 2.0)
pylab.grid()
pylab.axes().set_aspect('equal') # set the aspect ratio of the plot
fig = pylab.figure()
ax = fig.gca(projection='3d')
ax.set_aspect('equal')
pylab.plot(x_list_n3, y_list_n3, z_list_n3, '.')
pylab.title('Uniform sampling of the sphere surface')
ax.set_xlabel('$x$', fontsize=14)
ax.set_ylabel('$y$', fontsize=14)
ax.set_zlabel('$z$', fontsize=14)
pylab.show()
Populating the interactive namespace from numpy and matplotlib
Here, the renormalized distributions indicate that the gaussians are uniformly distibuted on the suface of the (hyper)sphere. We can also directly sample the points inside the 3d sphere of radius 1 as done by the following program.
%pylab inline
import random, math, pylab, mpl_toolkits.mplot3d
x_list, y_list, z_list = [],[],[]
nsamples = 1000
for sample in xrange(nsamples):
x, y, z = random.gauss(0.0, 1.0), random.gauss(0.0, 1.0), random.gauss(0.0, 1.0)
length = random.uniform(0.0, 1.0) ** (1.0 / 3.0) / math.sqrt(x ** 2 + y ** 2 + z ** 2)
x, y, z = x * length, y * length, z * length
#if z < 0.075 and z > -0.075 or z > 0.85 or z < -0.85:
x_list.append(x)
y_list.append(y)
z_list.append(z)
# graphics output
fig = pylab.figure()
ax = fig.gca(projection='3d')
ax.set_aspect('equal')
pylab.title('Uniform sampling inside the sphere\n(only shown three intervals for z)')
ax.set_xlabel('$x$', fontsize=14)
ax.set_ylabel('$y$', fontsize=14)
ax.set_zlabel('$z$', fontsize=14)
pylab.plot(x_list, y_list, z_list, '.')
pylab.show()
Populating the interactive namespace from numpy and matplotlib
The following program distributes the $d$ dimensional gaussian random variables on the surface of the $d$ dimensional hypersphere and calculates the corresponding radii distributions $\pi_d(r)$.
# %pylab inline
import random, math, pylab, numpy
nsamples = 100000
colors = ['0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9']
colors2 = ['r','g','b','y','c','m','0.7','0.8','0.9']
i = 0
for dimension in [2,4,6,10]:
r_samp = []
i += 1
#analytic distribution of r
dd = dimension
#plot the analytic distribution: a gamma distribution with k=d/2 and theta=2.
norm = math.gamma(dd / 2) * math.pow(2, (dd / 2)) / 2
x = [j * 0.1 for j in xrange(1, 400)]
gamma = [math.pow(r, (dd - 1)) * math.exp(- r ** 2 / 2) / norm for r in x]
#figure output
pylab.plot(x, gamma, color=colors2[i-1], label='analytic: d=%.0f' % dd)
#sampling distribution of r
for sample in range(nsamples):
R = [random.gauss(0.0, 1.0) for d in range(dimension)]
radius = math.sqrt(sum(x ** 2 for x in R))
r_samp.append(radius)
print numpy.var(r_samp) #calculate the variances for each dimension
#print [x / radius for x in R] #print the points that are rescaled to be on the surface of the sphere
pylab.hist(r_samp, bins=150, normed=True, color=colors[i], histtype='step', label='sampling: d=%.0f' % dimension)
pylab.legend()
pylab.xlim(0,7)
pylab.title('$r$ distributions in diffrent dimensions, $\sigma=1$')
pylab.xlabel('$r$', fontsize = 15)
pylab.ylabel('$\pi_d(r,\sigma=1)$', fontsize = 15)
pylab.show()
0.42951734083905774 0.46408913359880166 0.4790036671204991 0.4890500521524456
Using a similar strategy as we did for 1D, one can evaluate an $n$ dimensional gaussian integral by variable transformations and get $1=\left(\frac{1}{\sqrt{2\pi}}\right)^d\int_0^\infty \text{d}r r^{d-1} e^{-r^2/2} \int \text{d}\Omega_d$. Hence we see that the radii $r(d)$ in $d$ dimensions is distributed according to $\pi_d(r)\propto r^{d-1} e^{-r^2/2}$. In other words, the transformed variable square radius $\rho=r^2$ is distributed according to the gamma distirbution $\Gamma(\rho)\propto\rho^{\frac{d-2}{2}}e^{-\rho^2/2}$. This is confirmed by the following program. Furthermore, intriguingly we see that each distribution is approximately centered around the value $r^2=d$. That means that in order to renormalise a gaussian random sample in $d$ dimensions, it is most likely that a successful approach is to rescale it by about an amount $\sqrt{d}$. This implies that in very high dimensions, the renormalising points towards the surface of the unit hypersphare becomes unnecessary if the variance $\sigma$ of the Gaussian was choosen correctly.
# %pylab inline
import random, math, pylab
nsamples = 20000
colors = ['0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9']
colors2 = ['r','g','b','y','c','m','0.7','0.8','0.9']
i = 0
for dimension in [4,6,8,10]:
i += 1
dd = dimension
#plot the analytic distribution: a gamma distribution with k=d/2 and theta=2.
norm = math.gamma(dd / 2) * math.pow(2, (dd / 2))
x = [j * 0.1 for j in xrange(1, 400)]
gamma = [math.pow(rho, (dd / 2 - 1)) * math.exp(- rho / 2) / norm for rho in x]
#figure output
pylab.plot(x, gamma, color=colors2[i-1], label='analytic: d=%.0f' % dd)
r2 = []
for sample in range(nsamples):
R = [random.gauss(0.0, 1.0) for d in range(dimension)]
radius2 = sum(x ** 2 for x in R)
r2.append(radius2)
#print [x / radius for x in R] #print the points that are rescaled to be on the surface of the sphere
#figure output
pylab.hist(r2, bins=150, normed=True, color=colors[i], histtype='step', label='sampling: d=%.0f' % dimension)
pylab.legend()
pylab.title('$r^2$ distributions in diffrent dimensions with $\sigma=1$')
pylab.xlabel('$r^2$', fontsize = 15)
pylab.ylabel('$\pi_d(r^2,\sigma=1)=\Gamma_{k=d/2,\\theta=2}(r^2)$', fontsize = 15)
pylab.show()
In this case, we may rescale all points not by their actual distance from the origin, but instead by this central (most likely) value of $\sqrt{d}$, i.e. $x_i\rightarrow y_i=x_i\sqrt{d}$. This corresponds to transforming the gaussian distribution as $\pi(x_i)\propto e^{-x_i^2/2}\rightarrow \pi(y_i)\propto e^{-y_i^2/2\sigma^2}=e^{-x_i^2d/2}$, or in other words, taking the variance as $\sigma_x=1\rightarrow\sigma_y=\sqrt{\frac{1}{d}}$. Incredibly, as we increase $d$ the variance decreases, i.e. instead of true renormalisation, if we choose a fixed scaling factor of $\sqrt{d}$ the risk that we have the wrong scaling reduces with increasing number of dimensions!
# %pylab inline
import random, math, pylab, numpy
nsamples = 200000
colors = ['0.1','0.2','0.3','0.4','0.5','0.6','0.7','0.8','0.9']
colors2 = ['r','g','b','y','c','m','0.7','0.8','0.9']
i = 0
for dimension in [3,5,10,18]:
i += 1
#sampling
r2 = []
for sample in range(nsamples):
R = [random.gauss(0.0, math.pow(1.0 / dimension ,1.0/2.0)) for d in range(dimension)]
radius2 = (sum(x ** 2 for x in R))
r2.append(radius2)
print numpy.var(r2) #calculate the variances for each dimension
#print [x / radius for x in R] #print the points that are rescaled to be on the surface of the sphere
#figure output
pylab.hist(r2, bins=150, normed=True, color=colors[i], histtype='step', label='d=%.0f' % dimension)
pylab.legend()
pylab.xlim(0,5)
pylab.title('Rescaled $r$ distributions in diffrent dimensions, $\sigma=\sqrt{1/d}$')
pylab.xlabel('$r$', fontsize = 15)
pylab.ylabel('$\pi_d(r,\sigma=\sqrt{1/d})$', fontsize = 15)
pylab.show()
0.6645927747878635 0.40020680536587616 0.20079653557590216 0.11102109574157934
For a given kinetic energy $E_k$, the speeds of $N$ particles are constrained by the equation $E_k=\frac{m}{2}\sum_{i=0}^{N-1}v_i^2$, where e.g. in d dimensions, the speed is given in terms of the two components as $v_i^2=\sum_{j=1}^{d}(v_i^j)^2$, where $d$ is either $1$, $2$ or $3$. Therefore, the energy constraint implies that a legal set of velocities is a point on the surface of the $dN$ dimensional hypersphere of radius $r=\sqrt{\frac{2E_k}{m}}$. Hence the problem of finding the velocity distribution of a statistical ensemble is reduced to sampling a random point on the surface of a hypersphere.
On the other hand, since an ensemble contains about $6\times 10^{23}$ of particles, what we have is an hypersphere of a very large number of dimensions. Following the previous argument, for a sampled velocity component $v_i^j$, the suitable fixed scaling towards the surface of the hypersphere of radius $\sqrt{\frac{2E_k}{m}}$ is $\sqrt{\frac{2E_k}{mdN}}$ with certainty approaching to absolute as $N\rightarrow \infty$. (Please note that this implies that the kinetic energy is uniformly partitioned and distributed equally to each degree of freedom, i.e. to each velocity component of each particle.) In other words, instead of sampling velocities as gaussian variables each with variance $\sigma=1$, if we choose a variance that scales each velocity towards the surface of the kinetic energy hypersphere, then the probability that we are mistaken on the value a specific velocity component decreases by increasing $N$!
If we take the mean energy per degree of freedom at a given temperature $T$ as $\frac{E_k}{dN}=\frac{1}{2}k_BT$, where $k_B$ is the Boltzmann constant, we therefore arrive at the celebrated Maxwell velocity distribution $\pi(v_x)dv_x=\sqrt{\frac{m}{2\pi k_BT}}\exp{\left(-\frac{1}{2}\frac{mv_x^2}{k_B T}\right)}dv_x$. The speed distributions clearly have the gamma profiles similar to the ones shown previously. The equally famous Boltzmann distribution for energies $\pi(E)\propto\exp{\left(-\frac{E}{k_BT}\right)}$ is also readily obtained as a result of this approach.
A naive random number generator (NRNG): congruential linear random number generator.
Observe that $\text{idum}\in [0,m-1]$ because of the $\mod{m}$ operation. The method generates an erratic sequence which can contain all integers $\in [0,m-1]$. For a seed equal to e.g. 1000, the NRNG yields the same sequence of numbers at each run. However, if we set the initial seed "idum" to e.g. the system time, then the sequence of random numbers is different at every run.
Note that the output of the program is periodic, if you hit a number a certain time, you'll generate the same sequence as you did the first time. The program generates the same numbers it generated before.
m = 134456 #mode
n = 8121
k = 28411
idum = 1000 #seed (initial integer value)
for iteration in xrange(200000):
idum = (idum * n + k) % m #modulo operation (fold back to the range given by m) "congruential"
ran = idum / float(m) #divide by m to obtain a "random" number between 0 and 1.
#print idum, ran, iteration
Direct sampling of discrete and one-dimensional continuous systems is an important subject.
Here we sample 1D gaussian distribution through Markov-chain rejection sampling. (Calculating the area under the gaussian with the same method that is used in the calculation of $\pi$ in the first week.)
import random, math, pylab
x = 0.0 #start at x=0
delta = 0.5 #maximum step size
data = []
for k in range(50000):
x_new = x + random.uniform(-delta, delta) #propose a random move with step size between \pm delta
#if the randomly generated number between 0 and 1 is less than the acceptance rate, then accept the sample
if random.uniform(0.0, 1.0) < \
math.exp (- x_new ** 2 / 2.0) / math.exp (- x ** 2 / 2.0): #acceptance rate
x = x_new
data.append(x)
pylab.hist(data, 100, normed = 'True') #histogram of the sample
x = [a / 10.0 for a in range(-50, 51)]
y = [math.exp(- a ** 2 / 2.0) / math.sqrt(2.0 * math.pi) for a in x]
pylab.plot(x, y, c='red', linewidth=2.0)
pylab.title('Theoretical Gaussian distribution $\pi(x)$ and \
\nnormalized histogram for '+str(len(data))+' Markov-chain samples', fontsize = 18)
pylab.xlabel('$x$', fontsize = 30)
pylab.ylabel('$\pi(x)$', fontsize = 30)
pylab.savefig('plot_markov_gauss.png')
pylab.show()
Let us do the same sampling by direct sampling on a limited range given by a rectangle. (Note that we have to introduce an inelegant cutoff $\pm x_\text{cut}$ for the $x$ range of the rectangle.)
import random, math
y_max = 1.0 / math.sqrt(2.0 * math.pi)
x_cut = 5.0
n_data = 10000
n_accept = 0
data = []
while n_accept < n_data:
#select a random position within the rectangle
y = random.uniform(0.0, y_max)
x = random.uniform(-x_cut, x_cut)
if y < math.exp( - x **2 / 2.0)/math.sqrt(2.0 * math.pi): #check whether the sample is below the gaussian curve
n_accept += 1
data.append(x)
pylab.hist(data, 100, normed = 'True') #histogram of the sample
x = [a / 10.0 for a in range(-50, 51)]
y = [math.exp(- a ** 2 / 2.0) / math.sqrt(2.0 * math.pi) for a in x]
pylab.plot(x, y, c='red', linewidth=2.0)
pylab.title('Theoretical Gaussian distribution $\pi(x)$ and \
\nnormalized histogram for '+str(len(data))+' direct samples', fontsize = 18)
pylab.xlabel('$x$', fontsize = 30)
pylab.ylabel('$\pi(x)$', fontsize = 30)
pylab.savefig('plot_direct_gauss.png')
pylab.show()
Rejection sampling can be problematic It is not due to the presence of rejections per se, but because the rejection rate can become really enormous, even prohibitive. To see this, let us look at a function that is not quite as friendly as a Gaussian distribution: the distribution $\pi(x)=\frac{2}{\sqrt{x}}$. Here, we do not know what box size to choose: what should $x_\text{cut}$ be? With increasing $x_\text{cut}$, the rejection rate will increase enormously!
Below we compare the results from direct sampling with a cut-off and Markov-chain sampling methods.
import random, math, pylab
y_max = 100.0
x_cut = 1.0
n_data = 10000
data = []
n_accept = 0
while n_accept < n_data:
y = random.uniform(0.0, y_max)
x = random.uniform(0.0, x_cut)
if y < 1.0 / (2.0 * math.sqrt(x)):
n_accept += 1
data.append(x)
pylab.hist(data, bins=100, normed='True')
x = [a / 100.0 for a in xrange(1, 100)]
y = [1.0 / (2.0 * math.sqrt(a)) for a in x]
pylab.plot(x, y, 'red', linewidth = 2)
pylab.title('Theoretical distribution $\pi(x)={1}/{(2 \sqrt{x})}$ and normalized\
\n histogram for '+str(n_accept)+' accepted direct samples',fontsize=16)
pylab.xlabel('$x$', fontsize=18)
pylab.ylabel('$\pi(x)$', fontsize=18)
pylab.show()
import random, math, pylab
x = 0.2
delta = 0.5
data = []
y_max = 0
n_trials = 10000
for k in xrange(n_trials):
x_new = x + random.uniform(-delta, delta)
if x_new > 0.0 and x_new < 1.0:
if random.uniform(0.0, 1.0) < math.sqrt(x) / math.sqrt(x_new):
x = x_new
if 1.0 / math.sqrt(x) > y_max:
y_max = 1.0 / math.sqrt(x)
print y_max, x, k #print the maximum y values and their x values
data.append(x)
pylab.hist(data, bins=100, normed='True')
pylab.xlabel('$x$', fontsize=16)
pylab.ylabel('$\pi(x)$', fontsize=16)
x = [a / 100.0 for a in xrange(1, 101)]
y = [0.5 / math.sqrt(a) for a in x]
pylab.plot(x, y, linewidth=1.5, color='r')
pylab.title('Theoretical distribution $\pi(x)={1}/{(2 \sqrt{x})}$ and normalized\
\n histogram for '+str(len(data))+' Markov-chain samples',fontsize=16)
pylab.show()
2.2360679775 0.2 0 4.911152595 0.0414603671158 21 7.58973327104 0.0173598901546 68 11.3758766412 0.00772734410312 196 14.7582356921 0.00459125184976 369 29.4346418263 0.00115420377131 930 32.9788774849 0.000919450305395 3647 46.9006508302 0.000454613429938 6548 106.487244559 8.81870512329e-05 8023 122.781918606 6.63332336615e-05 8855
The tower sampling is rejection free but has logarithmic computational complexity and this is not optimal.
import random
# bisection search to find the bin corresponding to eta
def bisection_search(eta, w_cumulative):
kmin = 0
kmax = len(w_cumulative)
while True:
k = int((kmin + kmax) / 2)
if w_cumulative[k] < eta:
kmin = k
elif w_cumulative[k - 1] > eta:
kmax = k
else:
return k - 1
# sample an integer number according to weights
def tower_sample(weights):
sum_w = sum(weights)
w_cumulative = [0.0]
for l in xrange(len(weights)):
w_cumulative.append(w_cumulative[l] + weights[l])
eta = random.random() * sum_w
sampled_choice = bisection_search(eta, w_cumulative)
return sampled_choice
weights = [0.4, 0.3, 0.8, 0.1, 0.2]
n_samples = 20
for sample in xrange(n_samples):
print tower_sample(weights)s
Optimal.
import random, pylab
N = 5
pi = [[1.1 / 5.0, 0], [1.9 / 5.0, 1], [0.5 / 5.0, 2], [1.25 / 5.0, 3], [0.25 / 5.0, 4]]
x_val = [a[1] for a in pi]
y_val = [a[0] for a in pi]
pi_mean = sum(y_val) / float(N)
long_s = []
short_s = []
for p in pi:
if p[0] > pi_mean:
long_s.append(p)
else:
short_s.append(p)
table = []
for k in range(N - 1):
e_plus = long_s.pop()
e_minus = short_s.pop()
table.append((e_minus[0], e_minus[1], e_plus[1]))
e_plus[0] = e_plus[0] - (pi_mean - e_minus[0])
if e_plus[0] < pi_mean:
short_s.append(e_plus)
else:
long_s.append(e_plus)
if long_s != []:
table.append((long_s[0][0], long_s[0][1], long_s[0][1]))
else:
table.append((short_s[0][0], short_s[0][1], short_s[0][1]))
print table
samples = []
n_samples = 10000
for k in xrange(n_samples):
Upsilon = random.uniform(0.0, pi_mean)
i = random.randint(0, N-1)
if Upsilon < table[i][0]:
samples.append(table[i][1])
else: samples.append(table[i][2])
pylab.figure()
pylab.hist(samples, bins=N, range=(-0.5, N-0.5), normed=True)
pylab.plot(x_val, y_val,'ro', ms=8)
pylab.title("Histogram using Walker's method for a discrete distribution\n\
on $N=$"+str(N)+" choices ("+str(n_samples)+" samples)",fontsize=14)
pylab.xlabel('$k$',fontsize=20)
pylab.ylabel('$\pi(k)$',fontsize=20)
pylab.show()
[(0.05, 4, 3), (0.09999999999999998, 3, 1), (0.1, 2, 1), (0.17999999999999997, 1, 0), (0.19999999999999998, 0, 0)]
Let us investigate the continuum limit of tower sampling method.
For instance, for a gaussian distribution $\pi(x) \propto e^{-x^2}$, to obtain the position corresponding to a particular value $\Upsilon$ of the cumulative probability distribution $\Phi(x)\equiv \int^x_{-\infty}dx'\pi(x')$, we need to calculate the inverse of the error function $\text{erf}{(2\Phi-1)}$ which cannot be obtained analytically.
Below we pick a random value on the "tower" (i.e. we choose a value $\Upsilon$ of the cumulative distribution), and find the position that it corresponds to by numerically inverting the error function. The sampling is given as a histrogram.
import scipy.special, random, math
n_trials = 100000
data = []
for trial in xrange(n_trials):
Upsilon = random.uniform(0.0, 1.0) #pick a value for the cumulative distribution (or a value on the tower)
x = math.sqrt(2.0) * scipy.special.erfinv(2.0 * Upsilon - 1.0) #calculate x's from the numerical inverse of erf(x)
data.append(x)
pylab.hist(data, 100, normed = 'True') #histogram of the sample
x = [a / 10.0 for a in range(-50, 51)]
y = [math.exp(- a ** 2 / 2.0) / math.sqrt(2.0 * math.pi) for a in x]
pylab.plot(x, y, c='red', linewidth=2.0)
pylab.title('Theoretical Gaussian distribution $\pi(x)$ and \
\nnormalized histogram for '+str(len(data))+' "tower" samples', fontsize = 18)
pylab.xlabel('$x$', fontsize = 20)
pylab.ylabel('$\pi(x)\propto\exp{(-x^2/2)}$', fontsize = 20)
pylab.savefig('plot_tower_gauss.png')
pylab.show()
Let us now consider a divergent distribution $\pi(x)\propto x^\gamma$, $0<x\leq1$. For $\gamma < 0$, $\pi(x)$ blows up at $x=0$. However we have a finite cumulative distribution $\Phi(x)= x^{\gamma+1}$
import random
gamma = -0.5
n_trials = 10000
data = []
for trial in xrange(n_trials):
x = (random.uniform(0.0, 1.0)) ** (1.0 / (gamma + 1.0)) #simply invert the cumulative distribution
data.append(x)
pylab.hist(data, bins=100, normed='True')
pylab.xlabel('$x$', fontsize=16)
pylab.ylabel('$\pi(x)$', fontsize=16)
x = [a / 100.0 for a in xrange(1, 101)]
y = [0.5 / math.sqrt(a) for a in x]
pylab.plot(x, y, linewidth=1.5, color='r')
pylab.title('Theoretical distribution $\pi(x)={1}/{(2 \sqrt{x})}$ and normalized\
\n histogram for '+str(len(data))+' Markov-chain samples',fontsize=16)
pylab.show()
We found great rejection-free sampling programs for discrete and continuous one-dimensional distributions, and solved the direct sampling problem in these cases. Markov chain methods could generally be avoided. The situation is much more complicated in high-dimensional situations for example the hard-spheres simulations that we did last week or the quantum problem that we will turn our attention to next week.
We know a priori that the volume of the $d$-dimensional hypersphere ($d$-sphere) of radius $r$ is proportional to $r^{d}$, i.e. $V_d(r) \equiv V_\text{sph}(r,d)\propto r^{d}$. It follows that the surface area is given by $S_d(r)\equiv S_\text{sph}(r,d)=\frac{\text{d}V_\text{sph}(r,d)}{\text{d}r}=dC_nr^{d-1}$. Therefore, the task of calculating the volume is reduced to finding the proportionality constant $C_d$. If we consider $V_d(r\rightarrow \infty)$, by transforming the cartesian to the spherical coordinates $V_d(r\rightarrow \infty) = \int_0^{r \rightarrow \infty} \frac{\text{d}V_\text{sph}(d)}{\text{d}r'} \text{d}r' = C_d d \int_0^{r \rightarrow \infty} r'^{d-1} \text{d}r' = \int_{-\infty}^\infty \int_{-\infty}^\infty \cdots \int_{-\infty}^\infty \text{d}x_0\text{d}x_1 \cdots\text{d}x_{d-1} = \int \text{d}\Omega_d \int_{0}^\infty \text{d}r $, where $\text{d}\Omega_d$ is the $d$ dimensional solid angle differential, we get an integral relation given by $C_d d \int_0^{r \rightarrow \infty} r'^{d-1} \text{d}r' = \int \text{d}\Omega_d\int_{0}^\infty \text{d}r $.
Alternatively we consider the gaussian $e^{-r^2}$, where $r^2=x_0^2+\cdots+x_{d-1}^2$. We than have $I^d = \left( \int_{-\infty}^\infty e^{-x^2} \text{d}x \right)^d = \int \text{d}\Omega_d \int_{0}^\infty e^{-r^2} \text{d}r = C_d d \int_0^\infty r^{d-1}e^{-r^2} \text{d}r$ by the above equality. Hence it follows that $C_d=\frac{\left( \int_{-\infty}^\infty e^{-x^2} \text{d}x \right)^d}{d \int_0^\infty r^{d-1}e^{-r^2} \text{d}r}=\frac{2\pi^{d/2}}{\Gamma(d/2)}=\frac{2\pi^{d/2}}{\left(\frac{d}{2}-1\right)!}$. Using this, finally obtain the result that we were after: $S_d(r)=\frac{2\pi^{d/2}r^{d-1}}{\left(\frac{d}{2}-1\right)!}$ and $V_d(r)=\int_0^r S_d(r') \text{d}r' = \frac{\pi^{d/2}r^d}{\left(\frac{d}{2}\right)!}$. Therefore the volume of the unit $d$-sphere is $V_d(r=1) = V_\text{sph}(d) = \frac{\pi^{d/2}}{\left(\frac{d}{2}\right)!}$.
import math, pylab
from pylab import MaxNLocator
volume, dimension = [], []
def V_sph(dim):
return math.pi ** (dim / 2.0) / math.gamma(dim / 2.0 + 1.0)
for d in range(0,20):
dimension.append(d)
volume.append(V_sph(d))
pylab.plot(dimension, volume, 'bo-')
pylab.title('Volume of the unit hypersphere in $d$ dimensions')
pylab.xlabel('$d$', fontsize = 15)
pylab.ylabel('$V_{sph}(d)$', fontsize = 15)
pylab.xlim(0,20)
pylab.ylim(0,6)
for i in range(0,20):
pylab.annotate(round(volume[i],2), (dimension[i], volume[i]), xytext=(10, 0), ha='left',
textcoords='offset points')
pylab.grid()
pylab.show()
Although we are indeed able to calculate the exact expression, since the calculation of the gamma function $\Gamma(d)$ is difficult for large $d$, it is not entirely pointless to develop a numerical approach for this problem. We can numerically calculate $V_\text{sph}(d)$ from the ratio between the volumes of unit hyperspheres in successive dimensions $d$, $d+1$.
We first calculate $\langle Q(d=3) \rangle$ by a Markov-chain Monte Carlo program that samples points $(x, y, z)$ inside the unit cylinder and counts the samples that land inside the unit sphere, as well.
We see that the obtained value agrees with the analytical value of $Q(3)=\frac{4}{3}$ (see the table above).
import random
x, y, z = 0.0, 0.0, 0.0 #initial point
delta = 0.1
Q_avg_3 = 0.0 #initialise the average Q(3)
Q_3 = 4.0 / 3.0
n_trials = 4000
n_runs = 1000
for j in range(n_runs):
n_hits = 0
for i in range(n_trials):
del_x, del_y = random.uniform(-delta, delta), random.uniform(-delta, delta)
if (x + del_x) ** 2 + (y + del_y) ** 2 < 1.0: #sample points inside the unit disk
x, y = x + del_x, y + del_y
z = random.uniform(-1.0, 1.0) #sample the z coordinate (inside the unit cylinder)
#raise the number of counts if the sample is also inside the unit sphere
if x**2 + y**2 + z**2 < 1.0: n_hits += 1
Q_avg_3 += 2.0 * n_hits / float(n_trials) / n_runs #take the average
print '<Q(3)> =', Q_avg_3
print '<Q(3)>/Q(3)=<Q(3)>/(4/3) =', Q_avg_3 / Q_3 #compare with the actual value
<Q(3)> = 1.33341 <Q(3)>/Q(3)=<Q(3)>/(4/3) = 1.0000575
Now we do the same for $\langle Q(4) \rangle$ by a straightforward tweak:
import random, math
x, y, z, t = 0.0, 0.0, 0.0, 0.0 #initial point
delta = 0.1
Q_avg_4 = 0.0 #initialise the average Q(4)
Q_4= 3.0 * math.pi / 8.0
n_trials = 4000
n_runs = 1000
for j in range(n_runs):
n_hits = 0
for i in range(n_trials):
del_x, del_y, del_z = random.uniform(-delta, delta), random.uniform(-delta, delta), random.uniform(-delta, delta)
if (x + del_x)**2 + (y + del_y)**2 + (z + del_z)**2 < 1.0: #sample points inside the unit sphere
x, y, z = x + del_x, y + del_y, z + del_z
t = random.uniform(-1.0, 1.0) #sample the t coordinate (inside the unit cylinder)
#raise the number of counts if the sample is also inside the unit hypersphere
if x**2 + y**2 + z**2 + t**2 < 1.0: n_hits += 1
Q_avg_4 += 2.0 * n_hits / float(n_trials) / n_runs #take the average
print '<Q(4)> =', Q_avg_4
print '<Q(4)>/Q(4) =', Q_avg_4 / Q_4 #compare with the actual value
<Q(4)> = 1.175211 <Q(4)>/Q(4) = 0.997550079072
Obtain the volume of the 4-sphere $V_\text{sph}(4)$ from $\langle Q(4) \rangle$:
print 'Analytical: V_sph(4) = pi^2/2 =', math.pi**2 / 2 #analytical
print 'Numerical: V_sph(3)⋅<Q_4> =', 4.0 / 3.0 * math.pi * Q_avg_4
print 'Numerical: V_sph(2)⋅<Q_3>⋅<Q_4> = pi⋅<Q_3>⋅<Q_4> =', math.pi * Q_avg_3 * Q_avg_4
Analytical: V_sph(4) = pi^2/2 = 4.93480220054 Numerical: V_sph(3)⋅<Q_4> = 4.92681524536 Numerical: V_sph(2)⋅<Q_3>⋅<Q_4> = pi⋅<Q_3>⋅<Q_4> = 4.93103198336
Let us now efficiently generalise this program to sample uniformly distributed points inside the $d$-dimensional unit sphere, and calculate $\langle Q(d)\rangle$. Note that, here, instead of modifying all components of x at a time, as we did previously, modify only one component at each iteration $i$ (with $i=0, 1, 2,\cdots, n_\text{trials}$).
As seen previously, for small radii, e.g. $0<r<1$, the distribution of radii $\pi(r)$ is proportional to $r^{d-1}$ in random sampling inside the $d$-sphere. In the following, we verify our generalised program for $d=4$, $20$ by comparing the radius histograms with the analytic approximation to the radii distributions respectively given by $\pi_4(r) \approx 4 r ^3$ and $\pi_{20}(r) \approx 20 r ^{19}$ for small $r$.
%pylab inline
import pylab, random, math
#analytical formulae:
def V_sph(dim): #volume of "dim" dimensional sphere
return math.pi ** (dim / 2.0) / math.gamma(dim / 2.0 + 1.0)
def Q(dim):
return V_sph(dim) / V_sph(dim-1)
#monte carlo algorithm:
def markov_sphere(d, n_runs, n_trials):
x = [0.0] * d #initial point
delta = 0.1
Q_avg_d = 0 #initialise the average Q(d)
old_radius_square = sum(i**2 for i in x)
data = []
for j in range(n_runs):
n_hits = 0
if d == 1:
Q_avg_d = 1
break
for i in range(n_trials):
k = random.randint(0, d - 2)
x_old_k = x[k]
x_new_k = x_old_k + random.uniform(-delta, delta)
new_radius_square = old_radius_square + x_new_k ** 2 - x_old_k ** 2
if new_radius_square < 1.0: #check whether the position is inside the unit (d-1)-sphere
old_radius_square = new_radius_square #update the radius
x[k] = x_new_k #update a component
x[d-1] = random.uniform(-1, 1) #sample the d'th coordinate (inside the unit cylinder)
if old_radius_square + x[d-1]**2 < 1.0: #is x inside the unit d-sphere?
n_hits += 1 #raise the number of hits in case of accepted sample
data.append(math.sqrt(old_radius_square + x[d-1]**2)) #generate the radius histogram data
Q_avg_d += 2.0 * n_hits / float(n_trials) / n_runs #take the average of Q's
print 'Analytical: Q(d=%i)' % d, Q(d)
if d != 1:
print 'Numerical: <Q(d=%i)> =' % d, Q_avg_d
else:
print 'There is no cylinder in 1 dimension!'
return Q_avg_d, data
n_trials = 500
n_runs = 500
# 4-dimensional sphere radius distribution
Q_avg_d, data = markov_sphere(4, n_trials, n_trials)
pylab.hist(data, 100, normed = 'True') #histogram of the sample
x = [a / 100.0 for a in range(0, 100)]
y = [4 * a**3 for a in x]
pylab.plot(x, y, c='red', linewidth=2.0, label='$P_4(r) \sim 4 r^3, 0<r<1$')
pylab.xlabel('$x$', fontsize = 20)
pylab.ylabel('$\pi_4(x)$', fontsize = 20)
pylab.legend()
pylab.grid()
pylab.savefig('plot_radius_4.png')
pylab.show()
# 20-dimensional sphere radius distribution
Q_avg_d, data = markov_sphere(20, n_trials, n_trials)
pylab.hist(data, 100, normed = 'True') #histogram of the sample
x = [a / 100.0 for a in range(0, 100)]
y = [20 * a**19 for a in x]
pylab.plot(x, y, c='red', linewidth=2.0, label='$P_{20}(r) \sim 20 r^{19}, 0<r<1$')
pylab.xlabel('$x$', fontsize = 20)
pylab.ylabel('$\pi_{20}(x)$', fontsize = 20)
pylab.legend()
pylab.grid()
pylab.savefig('plot_radius_20.png')
pylab.show()
# 200-dimensional sphere radius distribution
Q_avg_d, data = markov_sphere(200, n_trials, n_trials)
pylab.hist(data, 100, normed = 'True') #histogram of the sample
x = [a / 100.0 for a in range(0, 100)]
y = [200 * a**199 for a in x]
pylab.plot(x, y, c='red', linewidth=2.0, label='$P_{200}(r) \sim 200 r^{199}, 0<r<1$')
pylab.xlabel('$x$', fontsize = 20)
pylab.ylabel('$\pi_{200}(x)$', fontsize = 20)
pylab.legend()
pylab.grid()
pylab.savefig('plot_radius_200.png')
pylab.show()
Populating the interactive namespace from numpy and matplotlib Analytical: Q(d=4) 1.1780972451 Numerical: <Q(d=4)> = 1.186448
Analytical: Q(d=20) 0.553539364154 Numerical: <Q(d=20)> = 0.551448
Analytical: Q(d=200) 0.177023967696 Numerical: <Q(d=200)> = 0.17836
We have now reached to the point where we can numerically calculate $V_\text{sph}(n)=Q(n) \cdots Q(3)Q(2) = 2 Q(n) \cdots Q(3)Q(2)$ for any dimension $n$ through successive calculations of $Q(d)$, $d\in [2,n]$.
import random, math
#Calculate the volume of the d_max dimensional sphere through the formula V(n)=Q(n)...Q(2)V(2)
def V_sph_markov(dim, n_runs, n_trials):
V = []
V_d = V_sph(1) #initialise the volume of "d_max" dimensional sphere
print 'The analytical value for the volume of the unit %i-sphere is =' % dim, V_sph(dim), '.'
print '___________________________________________'
for d in range(1, dim + 1):
Q_avg_d, data = markov_sphere(d, n_runs, n_trials)
V_d *= Q_avg_d
V.append(V_d) #save the volume for each d value along the way until d_max
print 'Analytical: V_sph(d=%i)' % d, V_sph(d)
print 'Numerical: V_sph(d=%i)' % d, V_d
print '___________________________________________'
print 'After %i runs' % n_runs, 'each consisting of %i trials' % n_trials, ', the numerical result for the volume of the unit %i-sphere (with analytical value' % dim, V_sph(dim), ') is found to be', V_d, '.'
return V
d_max = 200 #maximum dimension (dimension of the sphere in question)
n_trials = 1000 #each dimension iteration takes "n_trials" number of iterations
n_runs = 500
V = V_sph_markov(d_max, n_runs, n_trials)
The analytical value for the volume of the unit 200-sphere is = 5.55883284203e-109 . ___________________________________________ Analytical: Q(d=1) 2.0 There is no cylinder in 1 dimension! Analytical: V_sph(d=1) 2.0 Numerical: V_sph(d=1) 2.0 ___________________________________________ Analytical: Q(d=2) 1.57079632679 Numerical: <Q(d=2)> = 1.567504 Analytical: V_sph(d=2) 3.14159265359 Numerical: V_sph(d=2) 3.135008 ___________________________________________ Analytical: Q(d=3) 1.33333333333 Numerical: <Q(d=3)> = 1.33864 Analytical: V_sph(d=3) 4.18879020479 Numerical: V_sph(d=3) 4.19664710912 ___________________________________________ Analytical: Q(d=4) 1.1780972451 Numerical: <Q(d=4)> = 1.172344 Analytical: V_sph(d=4) 4.93480220054 Numerical: V_sph(d=4) 4.91991405849 ___________________________________________ Analytical: Q(d=5) 1.06666666667 Numerical: <Q(d=5)> = 1.058728 Analytical: V_sph(d=5) 5.26378901391 Numerical: V_sph(d=5) 5.20885077132 ___________________________________________ Analytical: Q(d=6) 0.981747704247 Numerical: <Q(d=6)> = 0.989828 Analytical: V_sph(d=6) 5.16771278005 Numerical: V_sph(d=6) 5.15586634128 ___________________________________________ Analytical: Q(d=7) 0.914285714286 Numerical: <Q(d=7)> = 0.913492 Analytical: V_sph(d=7) 4.72476597033 Numerical: V_sph(d=7) 4.70984265582 ___________________________________________ Analytical: Q(d=8) 0.859029241216 Numerical: <Q(d=8)> = 0.865424 Analytical: V_sph(d=8) 4.05871212642 Numerical: V_sph(d=8) 4.07601087057 ___________________________________________ Analytical: Q(d=9) 0.812698412698 Numerical: <Q(d=9)> = 0.808216 Analytical: V_sph(d=9) 3.29850890274 Numerical: V_sph(d=9) 3.29429720177 ___________________________________________ Analytical: Q(d=10) 0.773126317094 Numerical: <Q(d=10)> = 0.770212 Analytical: V_sph(d=10) 2.55016403988 Numerical: V_sph(d=10) 2.53730723637 ___________________________________________ Analytical: Q(d=11) 0.738816738817 Numerical: <Q(d=11)> = 0.725876 Analytical: V_sph(d=11) 1.88410387939 Numerical: V_sph(d=11) 1.84177042751 ___________________________________________ Analytical: Q(d=12) 0.708699124003 Numerical: <Q(d=12)> = 0.705664 Analytical: V_sph(d=12) 1.33526276885 Numerical: V_sph(d=12) 1.29967108696 ___________________________________________ Analytical: Q(d=13) 0.681984681985 Numerical: <Q(d=13)> = 0.672892 Analytical: V_sph(d=13) 0.910628754783 Numerical: V_sph(d=13) 0.874538277045 ___________________________________________ Analytical: Q(d=14) 0.658077758003 Numerical: <Q(d=14)> = 0.66044 Analytical: V_sph(d=14) 0.599264529321 Numerical: V_sph(d=14) 0.577580059691 ___________________________________________ Analytical: Q(d=15) 0.636519036519 Numerical: <Q(d=15)> = 0.627748 Analytical: V_sph(d=15) 0.381443280823 Numerical: V_sph(d=15) 0.362574727311 ___________________________________________ Analytical: Q(d=16) 0.616947898128 Numerical: <Q(d=16)> = 0.619088 Analytical: V_sph(d=16) 0.235330630359 Numerical: V_sph(d=16) 0.224465662782 ___________________________________________ Analytical: Q(d=17) 0.599076740253 Numerical: <Q(d=17)> = 0.603332 Analytical: V_sph(d=17) 0.140981106917 Numerical: V_sph(d=17) 0.135427317257 ___________________________________________ Analytical: Q(d=18) 0.582673014898 Numerical: <Q(d=18)> = 0.57796 Analytical: V_sph(d=18) 0.0821458866111 Numerical: V_sph(d=18) 0.0782715722821 ___________________________________________ Analytical: Q(d=19) 0.567546385503 Numerical: <Q(d=19)> = 0.563264 Analytical: V_sph(d=19) 0.0466216010301 Numerical: V_sph(d=19) 0.0440875588899 ___________________________________________ Analytical: Q(d=20) 0.553539364154 Numerical: <Q(d=20)> = 0.547916 Analytical: V_sph(d=20) 0.02580689139 Numerical: V_sph(d=20) 0.0241562789167 ___________________________________________ Analytical: Q(d=21) 0.540520367146 Numerical: <Q(d=21)> = 0.542652 Analytical: V_sph(d=21) 0.013949150409 Numerical: V_sph(d=21) 0.0131084530667 ___________________________________________ Analytical: Q(d=22) 0.528378483965 Numerical: <Q(d=22)> = 0.525592 Analytical: V_sph(d=22) 0.00737043094571 Numerical: V_sph(d=22) 0.00688969806424 ___________________________________________ Analytical: Q(d=23) 0.517019481618 Numerical: <Q(d=23)> = 0.51786 Analytical: V_sph(d=23) 0.00381065638685 Numerical: V_sph(d=23) 0.00356789903955 ___________________________________________ Analytical: Q(d=24) 0.5063627138 Numerical: <Q(d=24)> = 0.502836 Analytical: V_sph(d=24) 0.0019295743094 Numerical: V_sph(d=24) 0.00179406808145 ___________________________________________ Analytical: Q(d=25) 0.496338702353 Numerical: <Q(d=25)> = 0.495672 Analytical: V_sph(d=25) 0.000957722408823 Numerical: V_sph(d=25) 0.000889269314068 ___________________________________________ Analytical: Q(d=26) 0.486887224807 Numerical: <Q(d=26)> = 0.487064 Analytical: V_sph(d=26) 0.000466302805768 Numerical: V_sph(d=26) 0.000433131069187 ___________________________________________ Analytical: Q(d=27) 0.477955787451 Numerical: <Q(d=27)> = 0.475656 Analytical: V_sph(d=27) 0.000222872124721 Numerical: V_sph(d=27) 0.000206021391845 ___________________________________________ Analytical: Q(d=28) 0.46949839535 Numerical: <Q(d=28)> = 0.472048 Analytical: V_sph(d=28) 0.000104638104925 Numerical: V_sph(d=28) 9.72519859778e-05 ___________________________________________ Analytical: Q(d=29) 0.461474553401 Numerical: <Q(d=29)> = 0.463496 Analytical: V_sph(d=29) 4.82878227389e-05 Numerical: V_sph(d=29) 4.50759064928e-05 ___________________________________________ Analytical: Q(d=30) 0.453848448838 Numerical: <Q(d=30)> = 0.454848 Analytical: V_sph(d=30) 2.19153534478e-05 Numerical: V_sph(d=30) 2.05026859164e-05 ___________________________________________ Analytical: Q(d=31) 0.446588277485 Numerical: <Q(d=31)> = 0.442332 Analytical: V_sph(d=31) 9.78713994674e-06 Numerical: V_sph(d=31) 9.06899406678e-06 ___________________________________________ Analytical: Q(d=32) 0.439665684812 Numerical: <Q(d=32)> = 0.439904 Analytical: V_sph(d=32) 4.30306958703e-06 Numerical: V_sph(d=32) 3.98948676595e-06 ___________________________________________ Analytical: Q(d=33) 0.433055299379 Numerical: <Q(d=33)> = 0.43276 Analytical: V_sph(d=33) 1.86346708826e-06 Numerical: V_sph(d=33) 1.72649029283e-06 ___________________________________________ Analytical: Q(d=34) 0.426734341141 Numerical: <Q(d=34)> = 0.429012 Analytical: V_sph(d=34) 7.95205400148e-07 Numerical: V_sph(d=34) 7.4068505351e-07 ___________________________________________ Analytical: Q(d=35) 0.420682290826 Numerical: <Q(d=35)> = 0.422132 Analytical: V_sph(d=35) 3.34528829411e-07 Numerical: V_sph(d=35) 3.12666863008e-07 ___________________________________________ Analytical: Q(d=36) 0.414880609443 Numerical: <Q(d=36)> = 0.41544 Analytical: V_sph(d=36) 1.38789524622e-07 Numerical: V_sph(d=36) 1.29894321568e-07 ___________________________________________ Analytical: Q(d=37) 0.409312499182 Numerical: <Q(d=37)> = 0.416892 Analytical: V_sph(d=37) 5.68082871833e-08 Numerical: V_sph(d=37) 5.41519035072e-08 ___________________________________________ Analytical: Q(d=38) 0.403962698668 Numerical: <Q(d=38)> = 0.405004 Analytical: V_sph(d=38) 2.29484289973e-08 Numerical: V_sph(d=38) 2.1931737528e-08 ___________________________________________ Analytical: Q(d=39) 0.398817306895 Numerical: <Q(d=39)> = 0.397208 Analytical: V_sph(d=39) 9.15223065016e-09 Numerical: V_sph(d=39) 8.71146160003e-09 ___________________________________________ Analytical: Q(d=40) 0.393863631201 Numerical: <Q(d=40)> = 0.395044 Analytical: V_sph(d=40) 3.60473079746e-09 Numerical: V_sph(d=40) 3.44141063632e-09 ___________________________________________ Analytical: Q(d=41) 0.389090055507 Numerical: <Q(d=41)> = 0.3891 Analytical: V_sph(d=41) 1.40256490607e-09 Numerical: V_sph(d=41) 1.33905287859e-09 ___________________________________________ Analytical: Q(d=42) 0.384485925696 Numerical: <Q(d=42)> = 0.384944 Analytical: V_sph(d=42) 5.39266466261e-10 Numerical: V_sph(d=42) 5.15460371297e-10 ___________________________________________ Analytical: Q(d=43) 0.380041449565 Numerical: <Q(d=43)> = 0.383916 Analytical: V_sph(d=43) 2.0494360954e-10 Numerical: V_sph(d=43) 1.97893483907e-10 ___________________________________________ Analytical: Q(d=44) 0.375747609203 Numerical: <Q(d=44)> = 0.376504 Analytical: V_sph(d=44) 7.7007071306e-11 Numerical: V_sph(d=44) 7.45076882649e-11 ___________________________________________ Analytical: Q(d=45) 0.371596084019 Numerical: <Q(d=45)> = 0.373652 Analytical: V_sph(d=45) 2.86155261391e-11 Numerical: V_sph(d=45) 2.78399467355e-11 ___________________________________________ Analytical: Q(d=46) 0.367579182916 Numerical: <Q(d=46)> = 0.366528 Analytical: V_sph(d=46) 1.05184717169e-11 Numerical: V_sph(d=46) 1.02041199971e-11 ___________________________________________ Analytical: Q(d=47) 0.363689784359 Numerical: <Q(d=47)> = 0.361596 Analytical: V_sph(d=47) 3.82546071052e-12 Numerical: V_sph(d=47) 3.68976897447e-12 ___________________________________________ Analytical: Q(d=48) 0.359921283272 Numerical: <Q(d=48)> = 0.36048 Analytical: V_sph(d=48) 1.37686472804e-12 Numerical: V_sph(d=48) 1.33008791992e-12 ___________________________________________ Analytical: Q(d=49) 0.356267543862 Numerical: <Q(d=49)> = 0.353968 Analytical: V_sph(d=49) 4.90532214888e-13 Numerical: V_sph(d=49) 4.70808560837e-13 ___________________________________________ Analytical: Q(d=50) 0.352722857607 Numerical: <Q(d=50)> = 0.352916 Analytical: V_sph(d=50) 1.73021924584e-13 Numerical: V_sph(d=50) 1.66155874056e-13 ___________________________________________ Analytical: Q(d=51) 0.349281905747 Numerical: <Q(d=51)> = 0.350532 Analytical: V_sph(d=51) 6.04334275546e-14 Numerical: V_sph(d=51) 5.82429508447e-14 ___________________________________________ Analytical: Q(d=52) 0.34593972573 Numerical: <Q(d=52)> = 0.342352 Analytical: V_sph(d=52) 2.09063233531e-14 Numerical: V_sph(d=52) 1.99395907076e-14 ___________________________________________ Analytical: Q(d=53) 0.34269168111 Numerical: <Q(d=53)> = 0.343668 Analytical: V_sph(d=53) 7.16442309573e-15 Numerical: V_sph(d=53) 6.85259925929e-15 ___________________________________________ Analytical: Q(d=54) 0.339533434512 Numerical: <Q(d=54)> = 0.339988 Analytical: V_sph(d=54) 2.43256117999e-15 Numerical: V_sph(d=54) 2.32980151697e-15 ___________________________________________ Analytical: Q(d=55) 0.336460923272 Numerical: <Q(d=55)> = 0.337352 Analytical: V_sph(d=55) 8.18461780536e-16 Numerical: V_sph(d=55) 7.85963201352e-16 ___________________________________________ Analytical: Q(d=56) 0.333470337468 Numerical: <Q(d=56)> = 0.331608 Analytical: V_sph(d=56) 2.7293272616e-16 Numerical: V_sph(d=56) 2.60631685274e-16 ___________________________________________ Analytical: Q(d=57) 0.330558100057 Numerical: <Q(d=57)> = 0.332512 Analytical: V_sph(d=57) 9.02201234027e-17 Numerical: V_sph(d=57) 8.66631629338e-17 ___________________________________________ Analytical: Q(d=58) 0.327720848891 Numerical: <Q(d=58)> = 0.329124 Analytical: V_sph(d=58) 2.95670154285e-17 Numerical: V_sph(d=58) 2.85229268374e-17 ___________________________________________ Analytical: Q(d=59) 0.324955420395 Numerical: <Q(d=59)> = 0.327268 Analytical: V_sph(d=59) 9.6079619284e-18 Numerical: V_sph(d=59) 9.33464122023e-18 ___________________________________________ Analytical: Q(d=60) 0.322258834742 Numerical: <Q(d=60)> = 0.3227 Analytical: V_sph(d=60) 3.0962506153e-18 Numerical: V_sph(d=60) 3.01228872177e-18 ___________________________________________ Analytical: Q(d=61) 0.319628282356 Numerical: <Q(d=61)> = 0.32114 Analytical: V_sph(d=61) 9.8964926591e-19 Numerical: V_sph(d=61) 9.67366400109e-19 ___________________________________________ Analytical: Q(d=62) 0.317061111601 Numerical: <Q(d=62)> = 0.316496 Analytical: V_sph(d=62) 3.13779296345e-19 Numerical: V_sph(d=62) 3.06167596169e-19 ___________________________________________ Analytical: Q(d=63) 0.314554817556 Numerical: <Q(d=63)> = 0.315032 Analytical: V_sph(d=63) 9.87007893147e-20 Numerical: V_sph(d=63) 9.64525901563e-20 ___________________________________________ Analytical: Q(d=64) 0.312107031733 Numerical: <Q(d=64)> = 0.312644 Analytical: V_sph(d=64) 3.08052103827e-20 Numerical: V_sph(d=64) 3.01553235968e-20 ___________________________________________ Analytical: Q(d=65) 0.309715512671 Numerical: <Q(d=65)> = 0.312776 Analytical: V_sph(d=65) 9.5408515266e-21 Numerical: V_sph(d=65) 9.43186149332e-21 ___________________________________________ Analytical: Q(d=66) 0.307378137312 Numerical: <Q(d=66)> = 0.307188 Analytical: V_sph(d=66) 2.93264917062e-21 Numerical: V_sph(d=66) 2.89735466841e-21 ___________________________________________ Analytical: Q(d=67) 0.305092893079 Numerical: <Q(d=67)> = 0.305388 Analytical: V_sph(d=67) 8.9473041985e-22 Numerical: V_sph(d=67) 8.84817347476e-22 ___________________________________________ Analytical: Q(d=68) 0.302857870587 Numerical: <Q(d=68)> = 0.300768 Analytical: V_sph(d=68) 2.70976149705e-22 Numerical: V_sph(d=68) 2.66124743966e-22 ___________________________________________ Analytical: Q(d=69) 0.300671256947 Numerical: <Q(d=69)> = 0.301144 Analytical: V_sph(d=69) 8.14747395346e-23 Numerical: V_sph(d=69) 8.01418698968e-23 ___________________________________________ Analytical: Q(d=70) 0.298531329579 Numerical: <Q(d=70)> = 0.302312 Analytical: V_sph(d=70) 2.43227623203e-23 Numerical: V_sph(d=70) 2.42278489722e-23 ___________________________________________ Analytical: Q(d=71) 0.296436450511 Numerical: <Q(d=71)> = 0.298316 Analytical: V_sph(d=71) 7.21015332887e-24 Numerical: V_sph(d=71) 7.227554994e-24 ___________________________________________ Analytical: Q(d=72) 0.294385061112 Numerical: <Q(d=72)> = 0.296176 Analytical: V_sph(d=72) 2.12256142835e-24 Numerical: V_sph(d=72) 2.1406283279e-24 ___________________________________________ Analytical: Q(d=73) 0.292375677217 Numerical: <Q(d=73)> = 0.291912 Analytical: V_sph(d=73) 6.20585335048e-25 Numerical: V_sph(d=73) 6.24875096455e-25 ___________________________________________ Analytical: Q(d=74) 0.290406884611 Numerical: <Q(d=74)> = 0.292244 Analytical: V_sph(d=74) 1.80222253786e-25 Numerical: V_sph(d=74) 1.82615997688e-25 ___________________________________________ Analytical: Q(d=75) 0.288477334854 Numerical: <Q(d=75)> = 0.2891 Analytical: V_sph(d=75) 5.19900354536e-26 Numerical: V_sph(d=75) 5.27942849317e-26 ___________________________________________ Analytical: Q(d=76) 0.286585741392 Numerical: <Q(d=76)> = 0.287824 Analytical: V_sph(d=76) 1.48996028555e-26 Numerical: V_sph(d=76) 1.51954622662e-26 ___________________________________________ Analytical: Q(d=77) 0.284730875959 Numerical: <Q(d=77)> = 0.287028 Analytical: V_sph(d=77) 4.24237697249e-27 Numerical: V_sph(d=77) 4.36152314334e-27 ___________________________________________ Analytical: Q(d=78) 0.282911565221 Numerical: <Q(d=78)> = 0.282004 Analytical: V_sph(d=78) 1.20021750954e-27 Numerical: V_sph(d=78) 1.22996697251e-27 ___________________________________________ Analytical: Q(d=79) 0.281126687656 Numerical: <Q(d=79)> = 0.28234 Analytical: V_sph(d=79) 3.37413172925e-28 Numerical: V_sph(d=79) 3.4726887502e-28 ___________________________________________ Analytical: Q(d=80) 0.279375170655 Numerical: <Q(d=80)> = 0.283056 Analytical: V_sph(d=80) 9.42648627674e-29 Numerical: V_sph(d=80) 9.82965386876e-29 ___________________________________________ Analytical: Q(d=81) 0.277655987809 Numerical: <Q(d=81)> = 0.276744 Analytical: V_sph(d=81) 2.61732035873e-29 Numerical: V_sph(d=81) 2.72029773026e-29 ___________________________________________ Analytical: Q(d=82) 0.275968156379 Numerical: <Q(d=82)> = 0.278548 Analytical: V_sph(d=82) 7.22297074053e-30 Numerical: V_sph(d=82) 7.57733492167e-30 ___________________________________________ Analytical: Q(d=83) 0.274310734943 Numerical: <Q(d=83)> = 0.276316 Analytical: V_sph(d=83) 1.98133841231e-30 Numerical: V_sph(d=83) 2.09373887622e-30 ___________________________________________ Analytical: Q(d=84) 0.272682821184 Numerical: <Q(d=84)> = 0.272624 Analytical: V_sph(d=84) 5.40276947989e-31 Numerical: V_sph(d=84) 5.7080346739e-31 ___________________________________________ Analytical: Q(d=85) 0.271083549826 Numerical: <Q(d=85)> = 0.27008 Analytical: V_sph(d=85) 1.4646019295e-31 Numerical: V_sph(d=85) 1.54162600473e-31 ___________________________________________ Analytical: Q(d=86) 0.269512090705 Numerical: <Q(d=86)> = 0.269208 Analytical: V_sph(d=86) 3.94727928071e-32 Numerical: V_sph(d=86) 4.1501805348e-32 ___________________________________________ Analytical: Q(d=87) 0.267967646955 Numerical: <Q(d=87)> = 0.266368 Analytical: V_sph(d=87) 1.05774314073e-32 Numerical: V_sph(d=87) 1.10547528869e-32 ___________________________________________ Analytical: Q(d=88) 0.266449453311 Numerical: <Q(d=88)> = 0.269508 Analytical: V_sph(d=88) 2.8183508159e-33 Numerical: V_sph(d=88) 2.97934434106e-33 ___________________________________________ Analytical: Q(d=89) 0.264956774517 Numerical: <Q(d=89)> = 0.266176 Analytical: V_sph(d=89) 7.46741141638e-34 Numerical: V_sph(d=89) 7.93029959325e-34 ___________________________________________ Analytical: Q(d=90) 0.26348890383 Numerical: <Q(d=90)> = 0.264404 Analytical: V_sph(d=90) 1.96758004855e-34 Numerical: V_sph(d=90) 2.09680293365e-34 ___________________________________________ Analytical: Q(d=91) 0.26204516161 Numerical: <Q(d=91)> = 0.26082 Analytical: V_sph(d=91) 5.15594831803e-35 Numerical: V_sph(d=91) 5.46888141155e-35 ___________________________________________ Analytical: Q(d=92) 0.260624894005 Numerical: <Q(d=92)> = 0.263304 Analytical: V_sph(d=92) 1.34376848388e-35 Numerical: V_sph(d=92) 1.43997835119e-35 ___________________________________________ Analytical: Q(d=93) 0.259227471701 Numerical: <Q(d=93)> = 0.257688 Analytical: V_sph(d=93) 3.48341706628e-36 Numerical: V_sph(d=93) 3.71065141361e-36 ___________________________________________ Analytical: Q(d=94) 0.25785228875 Numerical: <Q(d=94)> = 0.25882 Analytical: V_sph(d=94) 8.98207063212e-37 Numerical: V_sph(d=94) 9.6039079887e-37 ___________________________________________ Analytical: Q(d=95) 0.256498761472 Numerical: <Q(d=95)> = 0.256448 Analytical: V_sph(d=95) 2.30388999259e-37 Numerical: V_sph(d=95) 2.46290299589e-37 ___________________________________________ Analytical: Q(d=96) 0.255166327409 Numerical: <Q(d=96)> = 0.257024 Analytical: V_sph(d=96) 5.87875148164e-38 Numerical: V_sph(d=96) 6.33025179615e-38 ___________________________________________ Analytical: Q(d=97) 0.253854444344 Numerical: <Q(d=97)> = 0.255 Analytical: V_sph(d=97) 1.49234719081e-38 Numerical: V_sph(d=97) 1.61421420802e-38 ___________________________________________ Analytical: Q(d=98) 0.252562589374 Numerical: <Q(d=98)> = 0.254532 Analytical: V_sph(d=98) 3.76911070755e-39 Numerical: V_sph(d=98) 4.10869170795e-39 ___________________________________________ Analytical: Q(d=99) 0.251290258037 Numerical: <Q(d=99)> = 0.251712 Analytical: V_sph(d=99) 9.47140802271e-40 Numerical: V_sph(d=99) 1.03420700719e-39 ___________________________________________ Analytical: Q(d=100) 0.25003696348 Numerical: <Q(d=100)> = 0.252564 Analytical: V_sph(d=100) 2.36820210188e-40 Numerical: V_sph(d=100) 2.61203458564e-40 ___________________________________________ Analytical: Q(d=101) 0.24880223568 Numerical: <Q(d=101)> = 0.249248 Analytical: V_sph(d=101) 5.89213977491e-41 Numerical: V_sph(d=101) 6.51044396403e-41 ___________________________________________ Analytical: Q(d=102) 0.247585620701 Numerical: <Q(d=102)> = 0.248128 Analytical: V_sph(d=102) 1.45880908343e-41 Numerical: V_sph(d=102) 1.61542343991e-41 ___________________________________________ Analytical: Q(d=103) 0.246386679994 Numerical: <Q(d=103)> = 0.246268 Analytical: V_sph(d=103) 3.59431126811e-42 Numerical: V_sph(d=103) 3.97827099699e-42 ___________________________________________ Analytical: Q(d=104) 0.245204989733 Numerical: <Q(d=104)> = 0.24584 Analytical: V_sph(d=104) 8.81343057595e-43 Numerical: V_sph(d=104) 9.78018141899e-43 ___________________________________________ Analytical: Q(d=105) 0.244040140185 Numerical: <Q(d=105)> = 0.243864 Analytical: V_sph(d=105) 2.15083083326e-43 Numerical: V_sph(d=105) 2.38503416156e-43 ___________________________________________ Analytical: Q(d=106) 0.242891735113 Numerical: <Q(d=106)> = 0.243004 Analytical: V_sph(d=106) 5.22419033025e-44 Numerical: V_sph(d=106) 5.79572841396e-44 ___________________________________________ Analytical: Q(d=107) 0.241759391211 Numerical: <Q(d=107)> = 0.242156 Analytical: V_sph(d=107) 1.26299707381e-44 Numerical: V_sph(d=107) 1.40347040981e-44 ___________________________________________ Analytical: Q(d=108) 0.240642737565 Numerical: <Q(d=108)> = 0.241412 Analytical: V_sph(d=108) 3.03931073379e-45 Numerical: V_sph(d=108) 3.38814598573e-45 ___________________________________________ Analytical: Q(d=109) 0.239541415145 Numerical: <Q(d=109)> = 0.241616 Analytical: V_sph(d=109) 7.28040794237e-46 Numerical: V_sph(d=109) 8.18630280489e-46 ___________________________________________ Analytical: Q(d=110) 0.238455076315 Numerical: <Q(d=110)> = 0.237576 Analytical: V_sph(d=110) 1.7360502315e-46 Numerical: V_sph(d=110) 1.94486907517e-46 ___________________________________________ Analytical: Q(d=111) 0.237383384378 Numerical: <Q(d=111)> = 0.239772 Analytical: V_sph(d=111) 4.12109479403e-47 Numerical: V_sph(d=111) 4.66325147893e-47 ___________________________________________ Analytical: Q(d=112) 0.236326013133 Numerical: <Q(d=112)> = 0.235712 Analytical: V_sph(d=112) 9.73921902418e-48 Numerical: V_sph(d=112) 1.0991843326e-47 ___________________________________________ Analytical: Q(d=113) 0.235282646463 Numerical: <Q(d=113)> = 0.235548 Analytical: V_sph(d=113) 2.29146922649e-48 Numerical: V_sph(d=113) 2.58910671175e-48 ___________________________________________ Analytical: Q(d=114) 0.234252977931 Numerical: <Q(d=114)> = 0.23726 Analytical: V_sph(d=114) 5.36783490142e-49 Numerical: V_sph(d=114) 6.14291458431e-49 ___________________________________________ Analytical: Q(d=115) 0.233236710407 Numerical: <Q(d=115)> = 0.231556 Analytical: V_sph(d=115) 1.25197615441e-49 Numerical: V_sph(d=115) 1.42242872948e-49 ___________________________________________ Analytical: Q(d=116) 0.232233555707 Numerical: <Q(d=116)> = 0.23398 Analytical: V_sph(d=116) 2.90750874e-50 Numerical: V_sph(d=116) 3.32819874125e-50 ___________________________________________ Analytical: Q(d=117) 0.231243234249 Numerical: <Q(d=117)> = 0.231212 Analytical: V_sph(d=117) 6.72341724645e-51 Numerical: V_sph(d=117) 7.69519487361e-51 ___________________________________________ Analytical: Q(d=118) 0.230265474726 Numerical: <Q(d=118)> = 0.23286 Analytical: V_sph(d=118) 1.54817086404e-51 Numerical: V_sph(d=118) 1.79190307827e-51 ___________________________________________ Analytical: Q(d=119) 0.229300013793 Numerical: <Q(d=119)> = 0.231596 Analytical: V_sph(d=119) 3.54995600478e-52 Numerical: V_sph(d=119) 4.14997585315e-52 ___________________________________________ Analytical: Q(d=120) 0.22834659577 Numerical: <Q(d=120)> = 0.229544 Analytical: V_sph(d=120) 8.10620368827e-53 Numerical: V_sph(d=120) 9.52602057235e-53 ___________________________________________ Analytical: Q(d=121) 0.227404972357 Numerical: <Q(d=121)> = 0.227776 Analytical: V_sph(d=121) 1.84339102565e-53 Numerical: V_sph(d=121) 2.16979886189e-53 ___________________________________________ Analytical: Q(d=122) 0.226474902362 Numerical: <Q(d=122)> = 0.225484 Analytical: V_sph(d=122) 4.1748180255e-54 Numerical: V_sph(d=122) 4.89254926574e-54 ___________________________________________ Analytical: Q(d=123) 0.225556151444 Numerical: <Q(d=123)> = 0.226272 Analytical: V_sph(d=123) 9.4165588681e-55 Numerical: V_sph(d=123) 1.10704690746e-54 ___________________________________________ Analytical: Q(d=124) 0.224648491859 Numerical: <Q(d=124)> = 0.226892 Analytical: V_sph(d=124) 2.11541574822e-55 Numerical: V_sph(d=124) 2.51180086927e-55 ___________________________________________ Analytical: Q(d=125) 0.223751702232 Numerical: <Q(d=125)> = 0.22336 Analytical: V_sph(d=125) 4.73327874594e-56 Numerical: V_sph(d=125) 5.6103584216e-56 ___________________________________________ Analytical: Q(d=126) 0.222865567321 Numerical: <Q(d=126)> = 0.222448 Analytical: V_sph(d=126) 1.054884853e-56 Numerical: V_sph(d=126) 1.24801301017e-56 ___________________________________________ Analytical: Q(d=127) 0.221989877805 Numerical: <Q(d=127)> = 0.222652 Analytical: V_sph(d=127) 2.34173759616e-57 Numerical: V_sph(d=127) 2.7787259274e-57 ___________________________________________ Analytical: Q(d=128) 0.221124430076 Numerical: <Q(d=128)> = 0.22286 Analytical: V_sph(d=128) 5.1781539134e-58 Numerical: V_sph(d=128) 6.1926686018e-58 ___________________________________________ Analytical: Q(d=129) 0.220269026039 Numerical: <Q(d=129)> = 0.22174 Analytical: V_sph(d=129) 1.14058691918e-58 Numerical: V_sph(d=129) 1.37316233576e-58 ___________________________________________ Analytical: Q(d=130) 0.219423472922 Numerical: <Q(d=130)> = 0.219764 Analytical: V_sph(d=130) 2.50271542977e-59 Numerical: V_sph(d=130) 3.01771647557e-59 ___________________________________________ Analytical: Q(d=131) 0.218587583092 Numerical: <Q(d=131)> = 0.219572 Analytical: V_sph(d=131) 5.4706251696e-60 Numerical: V_sph(d=131) 6.62606041973e-60 ___________________________________________ Analytical: Q(d=132) 0.217761173884 Numerical: <Q(d=132)> = 0.217624 Analytical: V_sph(d=132) 1.19128975882e-60 Numerical: V_sph(d=132) 1.44198977278e-60 ___________________________________________ Analytical: Q(d=133) 0.21694406743 Numerical: <Q(d=133)> = 0.217112 Analytical: V_sph(d=133) 2.58443245765e-61 Numerical: V_sph(d=133) 3.13073283549e-61 ___________________________________________ Analytical: Q(d=134) 0.216136090497 Numerical: <Q(d=134)> = 0.217476 Analytical: V_sph(d=134) 5.58589127551e-62 Numerical: V_sph(d=134) 6.8085925413e-62 ___________________________________________ Analytical: Q(d=135) 0.215337074338 Numerical: <Q(d=135)> = 0.218112 Analytical: V_sph(d=135) 1.20284948484e-62 Numerical: V_sph(d=135) 1.48503573637e-62 ___________________________________________ Analytical: Q(d=136) 0.214546854538 Numerical: <Q(d=136)> = 0.214568 Analytical: V_sph(d=136) 2.58067573454e-63 Numerical: V_sph(d=136) 3.18641147881e-63 ___________________________________________ Analytical: Q(d=137) 0.213765270876 Numerical: <Q(d=137)> = 0.21528 Analytical: V_sph(d=137) 5.51658847436e-64 Numerical: V_sph(d=137) 6.85970663158e-64 ___________________________________________ Analytical: Q(d=138) 0.212992167186 Numerical: <Q(d=138)> = 0.2138 Analytical: V_sph(d=138) 1.17499013463e-64 Numerical: V_sph(d=138) 1.46660527783e-64 ___________________________________________ Analytical: Q(d=139) 0.212227391229 Numerical: <Q(d=139)> = 0.211932 Analytical: V_sph(d=139) 2.49365090992e-65 Numerical: V_sph(d=139) 3.10820589742e-65 ___________________________________________ Analytical: Q(d=140) 0.211470794563 Numerical: <Q(d=140)> = 0.21104 Analytical: V_sph(d=140) 5.27334339284e-66 Numerical: V_sph(d=140) 6.55955772591e-66 ___________________________________________ Analytical: Q(d=141) 0.210722232426 Numerical: <Q(d=141)> = 0.213312 Analytical: V_sph(d=141) 1.11121069209e-66 Numerical: V_sph(d=141) 1.39923237763e-66 ___________________________________________ Analytical: Q(d=142) 0.209981563616 Numerical: <Q(d=142)> = 0.2126 Analytical: V_sph(d=142) 2.33333758631e-67 Numerical: V_sph(d=142) 2.97476803484e-67 ___________________________________________ Analytical: Q(d=143) 0.209248650381 Numerical: <Q(d=143)> = 0.209844 Analytical: V_sph(d=143) 4.88247740819e-68 Numerical: V_sph(d=143) 6.24237223503e-68 ___________________________________________ Analytical: Q(d=144) 0.208523358313 Numerical: <Q(d=144)> = 0.208628 Analytical: V_sph(d=144) 1.01811058604e-68 Numerical: V_sph(d=144) 1.30233363465e-68 ___________________________________________ Analytical: Q(d=145) 0.20780555624 Numerical: <Q(d=145)> = 0.209268 Analytical: V_sph(d=145) 2.11569036647e-69 Numerical: V_sph(d=145) 2.72536755056e-69 ___________________________________________ Analytical: Q(d=146) 0.207095116133 Numerical: <Q(d=146)> = 0.2099 Analytical: V_sph(d=146) 4.38149142144e-70 Numerical: V_sph(d=146) 5.72054648862e-70 ___________________________________________ Analytical: Q(d=147) 0.206391913001 Numerical: <Q(d=147)> = 0.207224 Analytical: V_sph(d=147) 9.04304396267e-71 Numerical: V_sph(d=147) 1.18543452556e-70 ___________________________________________ Analytical: Q(d=148) 0.205695824807 Numerical: <Q(d=148)> = 0.207192 Analytical: V_sph(d=148) 1.86011638667e-71 Numerical: V_sph(d=148) 2.45612550219e-71 ___________________________________________ Analytical: Q(d=149) 0.205006732377 Numerical: <Q(d=149)> = 0.205124 Analytical: V_sph(d=149) 3.81336382271e-72 Numerical: V_sph(d=149) 5.03810287512e-72 ___________________________________________ Analytical: Q(d=150) 0.204324519309 Numerical: <Q(d=150)> = 0.205616 Analytical: V_sph(d=150) 7.79163730025e-73 Numerical: V_sph(d=150) 1.03591456077e-72 ___________________________________________ Analytical: Q(d=151) 0.203649071897 Numerical: <Q(d=151)> = 0.203084 Analytical: V_sph(d=151) 1.58675970475e-73 Numerical: V_sph(d=151) 2.1037767266e-73 ___________________________________________ Analytical: Q(d=152) 0.20298027905 Numerical: <Q(d=152)> = 0.203152 Analytical: V_sph(d=152) 3.22080927656e-74 Numerical: V_sph(d=152) 4.27386449561e-74 ___________________________________________ Analytical: Q(d=153) 0.202318032212 Numerical: <Q(d=153)> = 0.203896 Analytical: V_sph(d=153) 6.51627794963e-75 Numerical: V_sph(d=153) 8.71423875198e-75 ___________________________________________ Analytical: Q(d=154) 0.20166222529 Numerical: <Q(d=154)> = 0.201592 Analytical: V_sph(d=154) 1.31408711193e-75 Numerical: V_sph(d=154) 1.75672081849e-75 ___________________________________________ Analytical: Q(d=155) 0.201012754584 Numerical: <Q(d=155)> = 0.20252 Analytical: V_sph(d=155) 2.64148270133e-76 Numerical: V_sph(d=155) 3.5577110016e-76 ___________________________________________ Analytical: Q(d=156) 0.200369518718 Numerical: <Q(d=156)> = 0.200548 Analytical: V_sph(d=156) 5.29272617567e-77 Numerical: V_sph(d=156) 7.13491825949e-77 ___________________________________________ Analytical: Q(d=157) 0.199732418568 Numerical: <Q(d=157)> = 0.200456 Analytical: V_sph(d=157) 1.05712899988e-77 Numerical: V_sph(d=157) 1.43023717463e-77 ___________________________________________ Analytical: Q(d=158) 0.199101357207 Numerical: <Q(d=158)> = 0.199004 Analytical: V_sph(d=158) 2.10475818619e-78 Numerical: V_sph(d=158) 2.84622918699e-78 ___________________________________________ Analytical: Q(d=159) 0.198476239835 Numerical: <Q(d=159)> = 0.200164 Analytical: V_sph(d=159) 4.17744490557e-79 Numerical: V_sph(d=159) 5.69712618985e-79 ___________________________________________ Analytical: Q(d=160) 0.197856973724 Numerical: <Q(d=160)> = 0.196636 Analytical: V_sph(d=160) 8.26536606916e-80 Numerical: V_sph(d=160) 1.12026010547e-79 ___________________________________________ Analytical: Q(d=161) 0.197243468159 Numerical: <Q(d=161)> = 0.197336 Analytical: V_sph(d=161) 1.63028946908e-80 Numerical: V_sph(d=161) 2.21067648172e-80 ___________________________________________ Analytical: Q(d=162) 0.19663563438 Numerical: <Q(d=162)> = 0.197848 Analytical: V_sph(d=162) 3.20573003977e-81 Numerical: V_sph(d=162) 4.37377920556e-81 ___________________________________________ Analytical: Q(d=163) 0.196033385532 Numerical: <Q(d=163)> = 0.197096 Analytical: V_sph(d=163) 6.28430112798e-82 Numerical: V_sph(d=163) 8.62054386299e-82 ___________________________________________ Analytical: Q(d=164) 0.19543663661 Numerical: <Q(d=164)> = 0.194248 Analytical: V_sph(d=164) 1.22818267589e-82 Numerical: V_sph(d=164) 1.6745234043e-82 ___________________________________________ Analytical: Q(d=165) 0.194845304408 Numerical: <Q(d=165)> = 0.195908 Analytical: V_sph(d=165) 2.39305627353e-83 Numerical: V_sph(d=165) 3.28052531089e-83 ___________________________________________ Analytical: Q(d=166) 0.194259307473 Numerical: <Q(d=166)> = 0.1949 Analytical: V_sph(d=166) 4.6487345444e-84 Numerical: V_sph(d=166) 6.39374383093e-84 ___________________________________________ Analytical: Q(d=167) 0.193678566058 Numerical: <Q(d=167)> = 0.195388 Analytical: V_sph(d=167) 9.00360240544e-85 Numerical: V_sph(d=167) 1.24926081964e-84 ___________________________________________ Analytical: Q(d=168) 0.193103002072 Numerical: <Q(d=168)> = 0.19258 Analytical: V_sph(d=168) 1.73862265395e-85 Numerical: V_sph(d=168) 2.40582648646e-85 ___________________________________________ Analytical: Q(d=169) 0.19253253904 Numerical: <Q(d=169)> = 0.192488 Analytical: V_sph(d=169) 3.34741433997e-86 Numerical: V_sph(d=169) 4.63092728726e-86 ___________________________________________ Analytical: Q(d=170) 0.19196710206 Numerical: <Q(d=170)> = 0.1929 Analytical: V_sph(d=170) 6.42593430237e-87 Numerical: V_sph(d=170) 8.93305873712e-87 ___________________________________________ Analytical: Q(d=171) 0.191406617759 Numerical: <Q(d=171)> = 0.191796 Analytical: V_sph(d=171) 1.22996635076e-87 Numerical: V_sph(d=171) 1.71332493354e-87 ___________________________________________ Analytical: Q(d=172) 0.190851014257 Numerical: <Q(d=172)> = 0.193368 Analytical: V_sph(d=172) 2.34740325544e-88 Numerical: V_sph(d=172) 3.3130221575e-88 ___________________________________________ Analytical: Q(d=173) 0.190300221124 Numerical: <Q(d=173)> = 0.191024 Analytical: V_sph(d=173) 4.46711358578e-89 Numerical: V_sph(d=173) 6.32866744613e-89 ___________________________________________ Analytical: Q(d=174) 0.189754169347 Numerical: <Q(d=174)> = 0.191204 Analytical: V_sph(d=174) 8.4765342785e-90 Numerical: V_sph(d=174) 1.21006653037e-89 ___________________________________________ Analytical: Q(d=175) 0.189212791289 Numerical: <Q(d=175)> = 0.191148 Analytical: V_sph(d=175) 1.6038687113e-90 Numerical: V_sph(d=175) 2.31301797147e-90 ___________________________________________ Analytical: Q(d=176) 0.188676020658 Numerical: <Q(d=176)> = 0.190292 Analytical: V_sph(d=176) 3.02611566105e-91 Numerical: V_sph(d=176) 4.40148815827e-91 ___________________________________________ Analytical: Q(d=177) 0.188143792469 Numerical: <Q(d=177)> = 0.1906 Analytical: V_sph(d=177) 5.69344876919e-92 Numerical: V_sph(d=177) 8.38923642967e-92 ___________________________________________ Analytical: Q(d=178) 0.187616043014 Numerical: <Q(d=178)> = 0.189652 Analytical: V_sph(d=178) 1.06818232918e-92 Numerical: V_sph(d=178) 1.59103546736e-92 ___________________________________________ Analytical: Q(d=179) 0.187092709829 Numerical: <Q(d=179)> = 0.187328 Analytical: V_sph(d=179) 1.99849126557e-93 Numerical: V_sph(d=179) 2.9804549203e-93 ___________________________________________ Analytical: Q(d=180) 0.186573731664 Numerical: <Q(d=180)> = 0.189 Analytical: V_sph(d=180) 3.72865973115e-94 Numerical: V_sph(d=180) 5.63305979936e-94 ___________________________________________ Analytical: Q(d=181) 0.186059048449 Numerical: <Q(d=181)> = 0.185216 Analytical: V_sph(d=181) 6.93750881567e-95 Numerical: V_sph(d=181) 1.0433328038e-94 ___________________________________________ Analytical: Q(d=182) 0.18554860127 Numerical: <Q(d=182)> = 0.187528 Analytical: V_sph(d=182) 1.28724505705e-95 Numerical: V_sph(d=182) 1.95654114031e-95 ___________________________________________ Analytical: Q(d=183) 0.185042332337 Numerical: <Q(d=183)> = 0.185428 Analytical: V_sph(d=183) 2.38194827645e-96 Numerical: V_sph(d=183) 3.62797510565e-96 ___________________________________________ Analytical: Q(d=184) 0.184540184959 Numerical: <Q(d=184)> = 0.184736 Analytical: V_sph(d=184) 4.39565175498e-97 Numerical: V_sph(d=184) 6.70217609117e-97 ___________________________________________ Analytical: Q(d=185) 0.184042103514 Numerical: <Q(d=185)> = 0.185384 Analytical: V_sph(d=185) 8.08984995301e-98 Numerical: V_sph(d=185) 1.24247621249e-97 ___________________________________________ Analytical: Q(d=186) 0.183548033427 Numerical: <Q(d=186)> = 0.184748 Analytical: V_sph(d=186) 1.48487604959e-98 Numerical: V_sph(d=186) 2.29544995304e-98 ___________________________________________ Analytical: Q(d=187) 0.183057921142 Numerical: <Q(d=187)> = 0.183552 Analytical: V_sph(d=187) 2.71818322792e-99 Numerical: V_sph(d=187) 4.21334429781e-99 ___________________________________________ Analytical: Q(d=188) 0.1825717141 Numerical: <Q(d=188)> = 0.183744 Analytical: V_sph(d=188) 4.96263371158e-100 Numerical: V_sph(d=188) 7.74176734657e-100 ___________________________________________ Analytical: Q(d=189) 0.182089360713 Numerical: <Q(d=189)> = 0.18248 Analytical: V_sph(d=189) 9.03642799993e-101 Numerical: V_sph(d=189) 1.4127177054e-100 ___________________________________________ Analytical: Q(d=190) 0.181610810341 Numerical: <Q(d=190)> = 0.18344 Analytical: V_sph(d=190) 1.64111301166e-101 Numerical: V_sph(d=190) 2.59148935879e-101 ___________________________________________ Analytical: Q(d=191) 0.181136013274 Numerical: <Q(d=191)> = 0.18094 Analytical: V_sph(d=191) 2.97264668265e-102 Numerical: V_sph(d=191) 4.68904084579e-102 ___________________________________________ Analytical: Q(d=192) 0.180664920704 Numerical: <Q(d=192)> = 0.182364 Analytical: V_sph(d=192) 5.37052977202e-103 Numerical: V_sph(d=192) 8.55112244802e-103 ___________________________________________ Analytical: Q(d=193) 0.180197484708 Numerical: <Q(d=193)> = 0.1805 Analytical: V_sph(d=193) 9.67755956468e-104 Numerical: V_sph(d=193) 1.54347760187e-103 ___________________________________________ Analytical: Q(d=194) 0.179733658226 Numerical: <Q(d=194)> = 0.180868 Analytical: V_sph(d=194) 1.73938318326e-104 Numerical: V_sph(d=194) 2.79165706895e-104 ___________________________________________ Analytical: Q(d=195) 0.179273395043 Numerical: <Q(d=195)> = 0.18082 Analytical: V_sph(d=195) 3.11825128544e-105 Numerical: V_sph(d=195) 5.04787431207e-105 ___________________________________________ Analytical: Q(d=196) 0.178816649766 Numerical: <Q(d=196)> = 0.177068 Analytical: V_sph(d=196) 5.57595247992e-106 Numerical: V_sph(d=196) 8.93817008689e-106 ___________________________________________ Analytical: Q(d=197) 0.178363377809 Numerical: <Q(d=197)> = 0.17626 Analytical: V_sph(d=197) 9.94545718822e-107 Numerical: V_sph(d=197) 1.57544185952e-106 ___________________________________________ Analytical: Q(d=198) 0.177913535373 Numerical: <Q(d=198)> = 0.18052 Analytical: V_sph(d=198) 1.76943144926e-107 Numerical: V_sph(d=198) 2.8439876448e-107 ___________________________________________ Analytical: Q(d=199) 0.177467079428 Numerical: <Q(d=199)> = 0.176668 Analytical: V_sph(d=199) 3.14015831549e-108 Numerical: V_sph(d=199) 5.02441609231e-108 ___________________________________________ Analytical: Q(d=200) 0.177023967696 Numerical: <Q(d=200)> = 0.176428 Analytical: V_sph(d=200) 5.55883284203e-109 Numerical: V_sph(d=200) 8.86447682334e-109 ___________________________________________ After 500 runs each consisting of 1000 trials , the numerical result for the volume of the unit 200-sphere (with analytical value 5.55883284203e-109 ) is found to be 8.86447682334e-109 .
As the following plot indicates, the Monte Carlo algorithm is able to calculate the volume of hyperspheres with exceptional accuracy up to very high dimensions.
import pylab
volume, dimension = [], []
for d in range(0,200):
dimension.append(d)
volume.append(V_sph(d))
pylab.semilogy(dimension, volume, 'b', label='Analytical', linewidth=4.0)
pylab.semilogy(dimension, V, 'r--', label='Monte Carlo', linewidth=3.0)
pylab.title('Volume of the unit hypersphere in $d$ dimensions')
pylab.xlabel('$d$', fontsize = 15)
pylab.ylabel('$V_{sph}(d)$', fontsize = 15)
pylab.ylim(0,10)
pylab.legend()
pylab.grid()
pylab.savefig('hypersphere_volumes.png')
pylab.show()
import sys, os
class HiddenPrints:
def __enter__(self):
self._original_stdout = sys.stdout
sys.stdout = open(os.devnull, 'w')
def __exit__(self, exc_type, exc_val, exc_tb):
sys.stdout.close()
sys.stdout = self._original_stdout
#Error calculation:
d = 20 #dimension of the sphere
n_runs = 10
Error = []
trials = [10, 100, 1000, 10000, 100000]
for trial in trials:
V_avg, V_avg_square = 0.0, 0.0
for run in range(n_runs):
with HiddenPrints():
v = V_sph_markov(d, 1, trial)
V_avg += v[d-1] / n_runs
V_avg_square += v[d-1] ** 2.0 / n_runs
print V_avg
Error.append(math.sqrt(V_avg_square - V_avg ** 2) / pow(n_runs, 0.5))
print Error
Dif = [] * len(trials)
for i in range(len(trials)):
Dif.append(abs(Error[i] - V_sph(d)))
print Dif
910059.1104 160100.804858 1.6837129243 0.0445860220085 0.0271548828196 [31358.716731116598, 13488.652683320564, 0.30053926397132036, 0.0030817988258178227, 0.0006134494292910471]